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Abstract. We propose a dedicated analysis approach for indirect Dark Matter searches with 
Imaging Air Cherenkov Telescopes. By using the full likelihood analysis, we take complete 
advantage of the distinct features expected in the gamma ray spectrum of Dark Matter origin, 
achieving better sensitivity with respect to the standard analysis chains. We describe the 
method and characterize its general performance. We also compare its sensitivity with that 
of the current standards for several Dark Matter annihilation models, obtaining gains of up 
to factors of order of 10. We compute the improved limits that can be reached using this 
new approach, taking as an example existing estimates for several benchmark models as well 
as the recent results from VERITAS on observations of the dwarf spheroidal galaxy Segue 
1. Furthermore, we estimate the sensitivity of Cherenkov Telescopes for monochromatic line 
signals. Predictions are made on improvement that can be achieved for MAGIC and CTA. 
Lastly, we discuss how this method can be applied in a global, sensitivity-optimized indirect 
Dark Matter search that combines the results of all Cherenkov observatories of the present 
generation. 
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1 Introduction 

The existence of Dark Matter (DM) has been confirmed by observational evidence on all 
scales, yet its nature still remains a mystery (for a review, see e.g. [1]). Among theories 
that try to describe DM and incorporate its presence in our image of the Universe today, 
the Cold Dark Matter paradigm offers the most satisfactory explanation. It requires the 
DM particle to be cold, neutral, stable on cosmological scales, consistent with the Big Bang 
nucleosynthesis and not excluded by the existing experimental constraints [2]. 

Among the possible candidates complying with these conditions, the Weakly Interacting 
Massive Particles (WIMPs) are the most widely considered. However, such particles do not 
exist within the framework of the Standard Model (SM), so one must go beyond its limits 
to look for WIMPs. For example, the Supersymmetric extension of the SM [3] suggests 
various natural DM particle candidates, with lightest neutralino being the most studied one. 
Assuming that neutralinos can self- annihilate, the resulting by-products are expected to be 
SM particles detectable from Earth, like electrons, positrons, photons and neutrinos. 

According to various models (see, e.g., [4]), the DM particle mass is expected to be in the 
few GeV - few TeV energy range; therefore, highly energetic photons resulting from WIMP 
annihilation might be visible to the Imaging Air Cherenkov Telescopes (IACTs). Gamma 
rays are especially attractive from the point of view of indirect DM searches: not only do 
they trace back to the place of their creation and can be detected from space and ground, 
their energy pattern also preserves information on the DM particle they originated from. 
Spectral features like cut-off [5], monochromatic line [6] or spectral hardening due to the 
internal bremsstrahlung [7] cannot be imitated by the conventional astrophysical sources; as 
such, they are considered to be the 'smoking guns' of DM detection. 

Cherenkov Telescopes have limited duty cycles and great variety of scientific objectives 
competing for the observation time. Their physics programs are primarily focused on detec- 
tion and study of astrophysical sources, with fundamental physics and cosmological issues 
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frequently left on the sidelines. As a consequence, standard analysis tools and methods used 
to process the observed data are usually optimized for sources with, in the majority of cases, 
featureless spectral distributions well described by a simple power law. 

Such analysis is suboptimal for DM searches, therefore, we propose an improved, ded- 
icated approach of optimized sensitivity for spectral features of DM origin. In comparison 
with some previous works, that addressed the issue by focusing only on the regions of the 
spectrum with the most peculiar, DM-induced features (see, e.g. [8] and references within), 
the method we are suggesting takes full profit of all spectral differences between the signal 
and the background. 

We describe the proposed approach in detail in section 2 and then go on to characterize 
it and compare its performance with respect to that of the standard IACT analysis in section 
3. In section 4, we apply the new method to several models of DM emission, show the 
improvement achievable with respect to the recently published experimental searches and 
make observability predictions for different annihilation channels and instruments. Lastly, 
section 5 is reserved for the discussion and concluding remarks. 



2 Full likelihood method 



In the standard analysis chain of IACTs, the existence of a source is established by a mere 
comparison of the integrated number of events detected in the source region (n) with the 
number of events detected from the control, background region(s) (m). Both n and m are 
random variables that obey Poisson statistics; therefore, the number of gamma-ray (g) and 
background (b) events present in the source region can be estimated by maximization of the 
following likelihood function [9]: 

n\ ml 

where r is the normalisation between the signal and background regions (e.g. ratio of their 
associated observation times). This, Poisson likelihood to which we also refer as the "conven- 
tional" likelihood approach, is what is currently used in the standard analysis of the IACT 
data. Whilst acceptable for sources of astrophysical origin, this method does not make any 
distinction of the potential features present in the gamma-ray spectrum, and as such, it is 
suboptimal for the DM searches. 

We propose the use of an alternative, more DM-oriented approach: by making an a 
priori assumption on the expected spectral shape (which is fixed and known for a given DM 
model), and including it in the maximum likelihood analysis, we can completely exploit the 
spectral information of the events from DM annihilation and achieve better sensitivity with 
respect to the conventional method. This full likelihood function has, for a given DM model 
M with parameters 9, the following form: 

Nf Nobs n obs 
C(N ES t,M(0)\Nobs,E 1 ,...,E Nobs )= gg-^ e~ N ^x JJ V(E i; M(0)), (2.2) 

2 = 1 

with Nobs(= n + m) and Nest denoting the total number of observed and estimated events, 
respectively, in source and background regions. 

V(Ei\ M{6)) is the value of the probability density function (PDF) of the event i with 
measured energy E^. In general, V may also depend on the measured arrival time and 
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direction of the photon, which would reflect on the sensitivity to gamma rays with distinct 
temporal and spatial structures, respectively. However, signals from DM annihilation are 
expected to be steady, allowing us to integrate out the time in our analysis. As for the 
spatial signatures, they may have a more important role for, e.g., galaxy clusters [10], as 
they can be predicted from halo simulations, although, usually with great uncertainties (see, 
for example, [11-13]). However, in this work we concentrate mainly on source-candidates 
that are of angular size smaller or comparable to the typical angular resolution of the IACTs 
(~ 0.1°) - hence, we do not expect a contribution from a likelihood function dependent on 
the direction, and we integrate it out as well. Therefore, we define the PDF as a function of 
measured energy only: 

HE; Mm= E P JSlBm , (2 . 3 ) 
j P(E-M(0))dE 

where E m i n and E max are the lower and upper limits of the considered energy range; 
P(E; M(6)) represents the differential rate of signal and background events, such that: 

P^M m = \ P ;f<\ ieB (2.4) 

with Pb(E) and Ps(E; M(0)) being the expected differential rates from the background (B) 
and source (S) regions, respectively: 

oo 

P B {E) = r J ^Rb(E; E')dE' (2.5) 



and 



P S (E; M(6)) = J ^R B (E; E>)dE> + J d l^M Rg{E , E < )dE > . (2 . 6) 


True energy is denoted with E'; d^s/dE' and d&c/dE' are the differential fluxes of cosmic 
(background) and gamma-ray (signal) radiations, and Rb(E; E') and Rq{E\E') are the 
telescope response functions to each of them. d&c/dE' contains the dependencies on the 
model parameters (6). 

However, in practice, Rb can be different for source and background regions, due to 
its dependence on the direction of the incoming particles within the observed field of view. 
Such discrepancies are measurable by the telescopes with relatively high precision, and the 
residual statistical and systematic uncertainties can be taken into account in the likelihood 
function through inclusion of the relevant nuisance parameters (see, e.g. [9, 14]). In this 
work we consider that Rb is equal in eq.(2.5) and eq.(2.6) and known with perfect precision. 
Nevertheless, in section 3.4 we evaluate the impact its uncertainties may have on our results. 

Apart from the shape of the signal spectral distribution, the given model M{6) also 
predicts the expected number of detected events for a given observation time Tobs'- 

Emax 

Nest = Tqbs J P(E; M(0))dE, (2.7) 
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Figure 1: Illustration of the advantage of the full likelihood method with respect to the conventional one. 
Red and orange lines show the assumed spectral energy distributions of the source and background regions, 
respectively, while the data points, with the same color code, represent the measured events (fine binning is 
used for demonstration purposes only, the full likelihood is unbinned). The levels of horizontal blue and cyan 
lines correspond to the average value within the energy range considered in the conventional method, with 
dots referring to the measurements. See the main text for more details. 



included in the full likelihood (eq.(2.2)) through the Poissonian term. 

Lastly, for the comparison of the full with conventional analysis, it is worthwhile to note 

that 

Emax 

b= Ioj3S J p B ( E ) dE (2.8) 



and 



1-smax 

g(0)=T OB s J Ps(E;M(6))dE-b. (2.9) 



Figure 1 illustrates the advantage of the full likelihood with respect to the conventional 
one. Both methods are based on comparisons of the collected data with the predictions 
from the signal and background models. The conventional method integrates the spectral 
information in a pre-optimized energy range (for details, see section 3.2), so that the only 
information used is that of the expected and measured number of events. On the other 
hand, the full likelihood compares the expected and measured energy distributions, thus fully 
profiting from the potential presence of DM spectral features. 



-4- 



3 Characterization of the full likelihood method 



In this section, in order to evaluate the performance of the full likelihood concept in the IACT 
analysis, we test it using fast simulations produced under a predefined set of conditions, and 
compare the results with those of the conventional method obtained under the exact same 
circumstances. 



3.1 The setup 

Response Function. The response functions Rb and Rg of an IACT are governed by 
its hardware design, reconstruction algorithms, selection criteria for quality of the events 
and for discrimination between the signal and background. They are computed by means 
of measurements and full Monte Carlo simulations, and typically can be represented as a 
product of three factors: effective area A e f /(£" ,p' ,t), angular (E(p;E',p',t)) and energy 
(G(E; E' ,p' ,t)) reconstruction functions, with p and p' referring to the measured and true 
directions of the incoming particle, and t to the time of the detection. As mentioned before, 
the spatial and temporal dependencies are integrated out in our analysis, so the response 
functions we use in this work have the form: 



Rb,g(E; E') = A effB G (E')Gb 7 g(E; E'). (3.1) 

The effective area A e ff is the area in which air showers can be observed by the instru- 
ment folded with the efficiency of all the cuts applied in the analysis. Cherenkov telescopes 
are not equally sensitive to gamma and cosmic ray showers, and effective areas for different 
particles are different as well. 

The energy resolution a is defined as the width of a Gaussian fit to the (E — E')/E' 
distribution, while the mean of that fit is the relative energy bias (i. However, in the real data 
analysis, the energy reconstruction function might need to be described by a more accurate 
parametrization. 

For the characterization of the full likelihood method in the following tests, as represen- 
tative response function of a current-generation IACT, we use the corresponding functions 
of the MAGIC Telescopes 1 [15]. 

Spectral Functions. The background emission is produced by the cosmic rays, with a 
flux well described by a simple power law: 

d<S> B „ w-« 



dE> 



AbE , (3.2) 



with spectral index a and intensity Ab- In practice, however, we only need the value of 
Pb(E) (eq.(2.5)), which is directly measured by the IACTs (or computed from Monte Carlo 
simulations for projected instruments). 

Regarding the spectral form of the signal emission, at this stage, we consider two simple 
cases: 

(a) power law (PL) of spectral index 7 and intensity Ap^: 

^ = A PL E'-y- (3.3) 



1 http://magic.mppmu.mpg.de 



- 5 - 





Figure 2: Contribution of the source region to the PDF before (purple) and after the convolution (green) with 
the response function of the telescope. Left: The spectral slope of a power law-shaped signal is harder after 
the convolution. Right: A monochromatic line is smoothed and widened due to the finite energy resolution. 
Shape of the background (left and right) is also affected by the response function. 



(b) a monochromatic line (L) at energy I and of intensity A^: 

d$ G 



dE< 



A L 5(E' - I). 



(3.4) 



The convolution of spectral and response functions yields the form of the PDF. As seen 
in figure 2, the original spectral shape is modified by the imperfect instrument, with features 
like line being smoothed and hardness of the power law being altered. 

Improvement Factor. In our tests, we choose the signal intensity (Api in the case of 
a PL, Al for the line-shaped signal) as a free parameter, whose value is to be estimated 
by maximization of the likelihood functions (eq.(2.1) and eq.(2.2)). Performance of the full 
likelihood with respect to the conventional one, for a given signal model M(0), is quantified 
by means of an Improvement Factor (IF) : 

IF{M(G)) = (CI cnm /CI full ), (3.5) 

i.e. the average ratio of the widths of the confidence intervals, CI cnvn and CIf u u, calculated 
by the corresponding methods, assuming a common confidence level. 

The Improvement Factor is, by construction, the improvement in the sensitivity of 
a given search expected by the use of the full likelihood over the conventional approach, 
provided that both methods produce unbiased estimators. We have explicitly checked this 
extreme, for several different models of signal emission, without finding any indications for 
the presence of bias (figure 3) . 

In this work, the confidence intervals are two-sided and computed following the HnC + 
1/2" rule and assuming one unconstrained degree of freedom. The maximization of the 
likelihood functions (eq.(2.1) and eq.(2.2)) is performed using the TMinuit class incorporated 
in the framework of ROOT [16, 17]. We have numerically confirmed that the obtained 
coverages are the expected ones. 

For the characterization of the method, the confidence intervals are calculated with 
95% confidence level, and their ratio averaged from 25 fast-simulated experiments. Each 
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Figure 3: Distribution of the free parameter values estimated by the conventional (blue) and full likelihood 
method (red), for a PL signal emission model, with 7 = 1.8. Test conditions are such that the expected 
parameter value is zero; results are obtained from 5000 fast-simulated experiments. 

simulation consists of 10 5 events 2 , randomly generated according to the PDF describing the 
expected background (i.e. the expected value of the signal intensity is zero), with r = 2, 
Emin = 100 GeV and E max = 10 TeV (unless specified otherwise). 

3.2 Optimization of the integration range 

By definition, the full likelihood takes complete advantage of the signal spectral information, 
so it makes sense to assume that maximal sensitivity with this method is achieved when the 
whole energy range is considered. For the conventional concept, however, this does not have 
to be the case, especially if some distinctive features are expected in the spectra. Here, we 
study the performance of each method for different energy integration ranges. For a chosen 
model and a given method, the optimal integration range is the one resulting in the best 
sensitivity. 

In the case of the power law-shaped signal, we fix one integration limit while varying 
the other: figure 4 shows the mean values of CIs, calculated with each method, for a signal 
of 7 = 1.8 spectral slope and the integration range of fixed E m i n (left) or fixed E max (right). 
As expected, in both cases, the full likelihood is best favoured when the entire energy range 
is considered. As for the conventional approach, the scenario with fixed E max and optimized 
Emin yields best sensitivity, and we shall always use such settings in the following tests. 

In the case of the spectral line, the sensitivity is optimized by restricting to those 
events in the vicinity of the peak. Figure 5 shows the CI widths of the full likelihood and 
conventional approaches, as a function of the integration range width (expressed in units 
of a), centered at I. Again, the full likelihood provides best constraints when the complete 

2 The number of events chosen for the characterization, for the selected setup, corresponds to ~200 hours of 
observations. This value depends very much on the chosen instrument and applied analysis cuts (in particular 
on the energy threshold), but does not have a significant role in the overall Improvement Factor value (for 
more details, see table 1, section 3.4) 
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Figure 4: Mean widths of the CIs, calculated with the conventional (blue) and full likelihood (red) methods, 
as a function of the integration range when E m in (left) or E max (right) is fixed. The considered signal model 
is a PL of spectral slope 7 = 1.8. Error bars are the RMS of the CI distributions. 
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Figure 5: Mean widths of the CIs, calculated with the conventional (blue) and full likelihood (red) methods, 
as a function of the integration range width given in units of a around the line energy I = 1 TeV. Error bars 
are the RMS of the CI distributions. 



energy span is integrated, while the conventional approach is most sensitive for a limited 
range. 

The Improvement Factor values given in the following sections are always calculated 
from the most constraining upper limits of both methods, using the whole energy range for 
the full likelihood and the optimized one for the conventional approach. 

3.3 Improvement Factor for different signal models 

In this section we compare the sensitivities of the full likelihood and conventional methods 
for various PL and line-shaped signals. 

Figure 6a shows the Improvement Factor as a function of the spectral slope 7 for PL 
models. In this example, for the case when 7 ~ 3.6, the shapes of signal and background 
differential rates are very alike, and therefore the improvement one gains from the use of the 
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Figure 6: Improvement Factor for different PL (a) and L (b) signal models (full line). Also shown are 
the optimal values of E m i n and integration range width (conventional approach) for the considered models 
(dashed line, right-hand axis). Error bars are the RMS of the IF distributions. 
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(a) (b) 

Figure 7: Improvement Factor as a function of spectral slope 7 (a) and of cutoff energy (b), for different signal 
models (full lines). Also shown are the optimal values of E m i n (conventional approach) for the considered 
models (dashed lines, right-hand axis) . Error bars are the RMS of the IF distributions. 



full likelihood is almost negligible. For harder spectral slopes the improvement on sensitivity 
of the full with respect to that of the conventional likelihood approach can reach up to 65%. 
The dashed line indicates the value of E m i n for which the conventional method yields the 
most constraining limit for the given model. For expected signal emissions of harder spectral 
indices, that dominate over the background radiation at higher energies, the conventional 
approach is optimized for the upper end of energy region. For increased 7, differences between 
signal and background concentrate at lower energies, so integration of the full energy range 
is preferred. 

For L models, depending on the line energy Z, the Improvement Factor can be between 
~40 and 65% (figure 6b). It is also interesting to note that the width of the optimal inte- 
gration range (in units of a at Z) for the conventional approach is almost constant for all the 
models and of order of 2.5 - 3. 

We can further elaborate the spectral shape of our signal by including additional features 
of physical interest. For example, the continuous, power law-shaped emission can abruptly 
cease at a certain energy, resulting in a sharp cutoff in the spectral distribution, smoothed 
by the response function of the detector. Figure 7a considers the case of PL models with 



- 9 - 



different spectral slopes 7 that all have a cutoff at a fixed energy of 1 TeV. In the presence 
of a cutoff, the Improvement Factor is lower than in the case of uninterrupted PL emission. 
This is especially noticeable for those signal models that dominate at high energies (7 < a), 
since their distinction from the background is partially erased by the cutoff. For the softer 
spectra, this effect is not that evident, as for those cases signal is more distinguishable from 
the background at lower energies, i.e. well below the cutoff. 

We also inspect how the Improvement Factor depends on the cutoff energy. Figure 7b 
shows that, for hard spectra, the higher the cutoff, the greater the improvement. For the 
soft spectra, Improvement Factor is enhanced by low-energy cutoffs to levels comparable to 
those obtained for spectral lines at similar energies. 

Lastly, we study the effect of adding a line to the PL-with-cutoff spectral distribution. 
For such models, we take the overall signal intensity as a free parameter, while the individual 
amplitudes of the PL, Api, and L, Al, are set in such a manner that the integrated areas 
corresponding to those emissions in the PDF are equal. As shown on figure 7a, presence of 
a line at the same energy as the cutoff (1 = 1 TeV) significantly boosts the Improvement 
Factor value, especially for soft spectra. Its contribution is obvious from the optimal E m in 
distribution as well: regardless of the value of 7, the most constraining limits from the 
conventional method are achieved when E m i n is just below the line, seeing how this feature 
is the one dominating the Improvement Factor value. 

3.4 Stability and robustness 

While in the previous subsection we varied the signal model, here we analyse the dependence 
of the Improvement Factor on the experimental parameters. Table 1 summarizes our results 
when, for several models of signal emission, we take different values of those parameters 
that are not affiliated with the DM model itself, but rather with the observational setup (r, 
number of events), energy resolution of the instrument (a) and energy range (E max ). For 
the majority of the considered settings, the observed variations of the Improvement Factor 
value are no more than 1% - 2%, with the following exceptions: 

• t: the greater the r, the lower the gain provided by the full likelihood method. This 
is especially noticeable for hard PL and L signal models; 

• a: in the case of the spectral line, the worse the a, the greater the improvement. It 
must be clarified, however, that this does not mean that a poor resolution yields more 
constraining upper limits, but that the advantage of the full likelihood approach is more 
significant; 

• E max : for the PL signal models of harder spectra, that dominate over the background 
at higher energies, the increase of energy range means also the greater Improvement 
Factor. On the other hand, for softer PL model, as well as for the considered L, change 
of this experimental parameter produces no significant effect. 

Additionally, we study the effect the presence of a signal in the data sample may have 
on the Improvement Factor value. Considering various signal intensities, that yield (for full 
likelihood) significances of up to 5 standard deviations, we find that the Improvement Factor 
increases up to ~10% for signal model of hard power law spectrum, whereas it has no sizable 
influence for other considered models (table 1). 

We also test the robustness of the full likelihood by assuming that the response function 
of the detector is not precisely known. For this, we simulate events with one response function, 
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Parameter 


Variation Range 

[units of the parameter] 


PL, 7 


= 1.8 


IF 

PL, 7 = 3.6 


L, I = 


1 TeV 


T- 


1-5 


1.91 - 


1.47 


1.02 - 1.01 


1.63 


- 1.26 


Number of events 


5 x 10 4 - 5 x 10 6 


1.66 - 


1.62 


1.03 - 1.02 


1.43 


- 1.41 


a [% of a magic] 


50 - 500 


1.65 - 


1.66 


1.01 - 1.11 


1.37 


- 2.83 


E ma x [TeV] 


10 - 50 


1.65 - 


1.82 


1.01 - 1.02 


1.40 


- 1.41 


Significance [std. dev] 


- 5 


1.65 - 


1.75 


1.01 - 1.01 


1.40 


- 1.42 



Table 1: Dependence of the Improvement Factor on different experimental parameters for three different 
representative signal models. 



Rq, but use a different one, R\y, f° r the likelihood maximization. Data are generated so that 
they contain a gamma-ray signal of intensity that yields a ha detection for Rq = Rw- We 
study how the significance of the detection by the full likelihood degrades when Rw ^ Ro ■ 

First, we consider the effect of using the wrong A e ff function, shifted for a fixed value in 
energy with respect to the real one. While for L signal models the sensitivity is not influenced 
by this discrepancy, for PL, especially those of soft spectral indices, the sensitivity decreases 
up to 5% for a 50 GeV shift. 

Next, we consider the case of unknown energy resolution a: for a power law-shaped 
signal, there is no significant effect - less than 1% decline in sensitivity for a factor 2 mistake 
in the estimate of a. On the other hand, in the case of a line, a a wrong by the same factor 
leads to a ~10% worse sensitivity. 

Lastly, we assume different energy bias fi functions for the simulations and for the 
likelihood analysis. Our findings show that, for [i values shifted from the actual ones by la 
at the given energy, the sensitivity of the analysis decreases ~5%. If the shift is 2a, the 
decline is ~20%. This means that, when searching for a line in the spectrum, we can take as 
up to la wide steps in our scan without risking a significant sensitivity degradation. 

Having in mind that even under these extreme and conservative conditions, the wors- 
ening in the sensitivity of the full likelihood is still smaller than the improvement one gains 
from its use (with respect to the conventional approach), we may conclude that this method 
is robust. 

As mentioned in section 2, the background in the source region may be known within 
some uncertainties. Here we estimate the effect this can have on the results of the conventional 
and full likelihood methods. 

First, we consider energy-dependent differences between the Rb functions in the source 
and background regions, parametrized as an extra power law term multiplying the first inte- 
gral in eq.(2.6). Its index is introduced in the likelihood functions as a nuisance parameter, 
with a Gaussian probability distribution of mean and width 0.04 (so that maximum de- 
viation of 5% is achieved at any energy). This results in the sensitivity decrease for both 
the full likelihood and conventional method, but more drastically for the latter one: for the 
case of the L models as well as the hard power law-shaped spectra, results from conventional 
approach are up to ~50% less constraining. For the full likelihood, the corresponding sen- 
sitivity losses are smaller: ~5% for L and ~25% for the PL signal models. Soft power law 
spectra are not affected (less than 1%), for either of the analysis methods. 
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The case of global (normalization) differences is examined by treating r as a nuisance 
parameter, with a Gaussian probability distribution of a 5% width. This leads to significant 
sensitivity losses for the conventional method: ~30% for L models and ~10% for hard PL 
signals. The full likelihood is again far more robust, exhibiting almost negligible worsening 
- less than 2% for both kinds of signal models. On the other hand, soft PL models result 
problematic for both methods, especially when the spectral shape of the signal is similar to 
that of the background. The conventional approach suffers from up to a factor ~8 worse 
sensitivity, also for all softer signal models. In the case of the full likelihood this is less 
pronounced (up to a factor ~4 sensitivity worsening), and its power is recovered as soon 
as the shape of the expected signal becomes different from that of the background. This is 
caused, for both methods, by high correlation (up to 0.99) between r and signal intensity 
when the signal and background are of similar spectral shapes. For other signal models the 
correlation is low, due to the energy range optimization applied in the conventional approach 
and the presence of the spectral term in the full likelihood. 



4 Sensitivity of the full likelihood method for Dark Matter searches 

So far we have characterized the performance of the full likelihood in a rather general way, 
by assuming generic spectral shapes. In this section, however, we explore its sensitivity for 
specific DM models, for the observations with MAGIC and CTA 3 . The following cases are 
considered: 

a) benchmark models (BM), as defined by Battaglia et al. (2009) [18] and by Bringmann, 
Doro and Fornasa (2009) [19]. We compare our results to those presented in [19]; 

b) secondary gamma rays from annihilation into SM particles (in particular, we study the 
bb, t + t~ and W + W~ channels). We also show the sensitivity improvement obtainable 
through the use of the full likelihood with respect to the recently published VERITAS 
results of Segue 1 observations [20] ; 

c) annihilation into 77. We compare the expected sensitivities to the possible hint of 
monochromatic line found in Fermi data by Weniger (2012) [21] and more recently by 
Su and Finkbeiner (2012) [22]. 

To compute the sensitivity of the full likelihood, we use, in eq.(2.2), the differential 
gamma-ray flux from the annihilation of DM particles, given as a product of two terms: 



dE< dE 

p 1 

■G 

,PP 



The particle physics term, d&Q P /dE', describes the characteristics of the chosen DM model: 



dW _ A. M dNa (4 2) 



dE' 4tt 1m\ dE' 



where m x refers to the DM particle mass, (av) is the thermally averaged cross-section, and 
dNc/dE' represents the differential gamma-ray rate per annihilation, summing all possible 
final states weighted by their corresponding branching ratios. 



https://www.cta-observatory.org 



- 12 - 



The effective astrophysical factor, J(Af2), depends on the distance and morphology of 
the source, and it is defined as the integral along the line of sight (los) of the squared DM 
density p, integrated over the solid angle AO of the signal region: 

J (AD;) = I dn f p 2 (r)ds. (4.3) 
J An Jios 

It is reasonable to assume that the J factor, for a given source and assumed DM distri- 
bution profile, is the same for every IACT of the current generation, since their point-spread 
functions are very alike, and therefore the signal region spans over similar solid angles AO. 
Unlike J, d&Q P jdE' does not depend on the observed source - its value is completely deter- 
mined for a given theoretical framework. 

In the upcoming subsections, we express the sensitivity of the full likelihood approach 
as the value of (crv) (taken as free parameter in the maximization of the likelihood) for which 
we would obtain a detection with a given statistical significance in a given Tobs- 

4.1 Benchmark models 

Bringmann, Doro and Fornasa (2009) [19] made observability predictions (requiring a 5cr 
detection in 50 hours) for two dwarf Spheroidal galaxies, Draco and Willman 1, for the case 
of several mSUGRA [23] BM models, and observations with MAGIC and CTA (although, the 
response functions attributed to each instrument are rather simplified and slightly optimistic) . 
In their calculations, they relied on the conventional likelihood approach, and made two 
studies: one, for which E m i n is the actual energy threshold of the analysis (70 GeV for 
MAGIC, 30 GeV for CTA), and the other, for which E m i n is optimized for each model based 
on the sensitivity curves of the instruments. In both cases, E max is selected as the DM 
particle mass m x . 

We have computed the sensitivities of the full and conventional likelihood approaches 
under the same circumstances studied in [19]: considering the same DM candidate sources, 
same BM models of DM emission, and same observatories (but with more realistic response 
functions: the actual one of MAGIC [15] and one of the latest estimates of the response func- 
tion of the CTA [24]). The results are shown in table 2, together with the basic characteristics 
of each of the studied BM models. Improvement Factors IF\ and IF2 represent the gain the 
full likelihood provides over the conventional method, for the two cases of integration ranges 
considered in [19]. 

The lowest Improvement Factors (although of values higher than 25%) are obtained, 
for both MAGIC and CTA, for the practically featureless, soft spectra of the model K' , as 
well as for the model /' from the bulk region, that has a cutoff at low energies. On the other 
hand, the greatest improvements are achieved in the case of the model BM4, characterized 
by the massive DM particle and hard spectrum. Models from coannihilation region, with 
particularly large internal bremsstrahlung contributions, J' and BM3, also show significant 
gain from the use of the full likelihood (above 60%). Despite these high improvements, 
however, estimated (crv) limits are still ~4 and ~3 orders of magnitude, for MAGIC and 
CTA respectively, away from the predicted values of these BM models. 

The fact that our results for (crv) f u ii are about factor ~2 less constraining than the 
conventional limits presented in [19], can be understood by taking into account that the latter 
were obtained assuming a somewhat idealized situation, with perfectly known background 
and with flat, optimistic response functions, while we consider circumstances of the real 
experiment and the actual (or latest from the simulations) responses of the detectors. 
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MAGIC (>70 GcV) CTA (>30 GeV) 



BM 


m x 
[GeV] 


[cm 3 s~ 


l ] 


[cm 3 s _1 


] 


m 


IF 2 


[cm 3 s _1 


] 


m 


IF 2 


r 


141 


3.6 x 10" 


-27 


5.65 x 10" 


23 


1.62 


1.57 


1.39 x 10" 


23 


1.48 


1.48 


./' 


316 


3.2 x 10" 


-28 


1.01 x 10" 


2:', 


3.64 


1.80 


1.91 x 10- 


24 


5.18 


1.65 


K' 


565 


2.6 x 10" 


-26 


3.91 x 10" 


2--! 


1.23 


1.23 


8.39 x 10- 


24 


1.58 


1.58 


BM3 


233 


9.2 x 10" 


-29 


7.21 x 10" 


25 


4.14 


1.89 


1.35 x 10" 


25 


6.62 


1.61 


BM4: 


1926 


2.6 x 10" 


-27 


2.87 x 10" 


23 


2.10 


2.10 


4.82 x 10- 


24 


3.81 


3.81 



Table 2: Characteristics of the studied BM models (mass m x and predicted annihilation cross section today 
av\ v= o), together with the upper limits on the (ov) value calculated with full likelihood method ((crv) f u u), 
for Willman 1 observations with MAGIC and CTA. We also quote the Improvement Factors obtainable from 
full likelihood with respect to the conventional approach, computed according to the prescription presented 
in [19]: IF\ is calculated for an integration range (for the conventional method) from energy threshold to m x , 
while for IF2 the integration is done from optimized lower limit to m x . 

4.2 Secondary gamma rays from annihilation into SM particles 

One of the most recent results from indirect DM searches with IACTs comes from the VERI- 
TAS Collaboration 4 , from observations of the dwarf spheroidal galaxy Segue 1 [20]. Although 
this source is considered to be one of the most DM dominated objects known so far [25, 26], 
nearly 50 hours of data showed no significant gamma-ray excess. Consequently, upper limits 
to the (crv) were derived, for the full energy range and relaying on the conventional likelihood 
analysis. 

We have computed how these results could be further strengthened if the full likelihood 
was used instead. For three final state channels, bb, t + t~ and W + W~ (100% branching 
ratios), we estimate the limits using both methods and calculate the Improvement Factor 
values for the upper limits on the DM particle annihilation cross section value (crv) . Following 
the prescription from [20], we assume r = 1/0.084, Tobs = 47.8 hours, J = 7.7xl0 18 GeV 2 
cm -5 , and we calculate the upper limits with 95% confidence level. In the lack of the 
VERITAS response function used for the analysis in [20], we use the A e ff function from 
Wood (2010) [27], a from [28] and we assume a conservative \i = 0. The total background rate 
is taken from [20], approximating the shape of Pg dependence on energy by that of MAGIC. 
For the dNc /dE' of the studied channels we use the parametrization from Cembranos et al. 
(2010) [29]. 

Figure 8a shows our estimates of (crv) limits, calculated by both full and conventional 
likelihood approaches, assuming observations with VERITAS, for considered channels and 
following the analysis prescription as given in [20]. The more massive the DM particle, 
the greater the improvement, especially for the r + r~ channel whose spectra gets harder for 
higher m x values. This is partially due to the non-optimization of the integration range in the 
VERITAS analysis. Gain achievable through the use of the full likelihood is quite significant: 
for example, the Improvement Factors for m x = 100 GeV, 1 TeV and 10 TeV, are 1.2, 1.6 and 
3.4 (for the bb channel), 1.2, 2.9 and 10.1 (for the r+r") and 1.3, 1.5 and 4.5 (for the W + W~ 
channel), respectively. Our estimates do not rely on the actual response function used in 
[20], and our limits from the conventional likelihood are slightly less constraining than those 
reported by VERITAS. Nevertheless, we confirm that the consistent Improvement Factor 

4 http://veritas. sao.arizona.edu 
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m z [GeV] m z [GeV] 

(c) (d) 

Figure 8: 95% confidence level upper limits on the {av) as a function of m x , considering different final 
state particles, and ~50 h of observations of the Segue 1 galaxy, (a) Results for VERITAS, from this work, 
estimated by the conventional (light green) and from the full likelihood approach (dark green). Exclusion 
lines for (b) X X -> bb, (c) X X t + t~ and (d) X X -> W + W~ channels, for MAGIC (blue) and CTA (red), 
obtained from both the conventional (dashed line) and full likelihood approach (full line), (b-d) Results from 
VERITAS, as given in [20], are also plotted for the comparison purposes (green). 



results (< 5% difference) are obtained when different A e ft function is used (McCutcheon 
(2012) [30]). From this, we infer the validity of the obtained Improvement Factor values also 
for the response function actually applied in the analysis from [20]. 

Additionally, we study the sensitivities of MAGIC and CTA for the gamma-ray spectra 
from bb, t + t~ and W + W~ channels, assuming the same observational and analysis condi- 
tions as in the case of VERITAS, but using the actual/simulated response functions of these 
instruments. Exclusion lines, calculated by means of both full likelihood and conventional 
methods, are shown on figure 8b-8d. As expected, the CTA results are always better than 
those of MAGIC: at lower energies, by a factor ~5, and by more than one order of magnitude 
at high energies. Again, the constraints from the full likelihood approach are more signifi- 
cantly improved with respect to the conventional ones for more massive DM particles, and 
the lowest (av) limits are achieved for the t + t~ channel (hardest spectrum). For MAGIC, 
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Improvement Factors for DM particle masses of m x = 100 GeV, 1 TeV and 10 TeV are rather 
relevant: 1.2, 1.4 and 3.4 (66), 1.3, 2.9 and 15.8 (t + t~) and 1.3, 1.5 and 4.5_ (W + W~). In 
the case of CTA, these values are even more significant: 1.3, 1.7 and 6.0 (66 channel), 1.3, 
4.8 and 33.1 (t + t~ channel) and 1.4, 2.0 and 8.5 (Ty+VF - channel), respectively. 

According to the results shown on figure 8b-8d, the sensitivity gain of the CTA with 
respect to VERITAS would be marginal, or even nonexistent, for certain annihilation channels 
and mass ranges. We have traced this inconsistency down to a probable overestimation of 
the VERITAS performance assumed in [20]. For that, we have used the response functions 
assumed for VERITAS and those of MAGIC and CTA to compute the integral sensitivity (5<r 
significance in 50 hours of observations) for a Crab-like spectrum 5 at the analysis threshold 6 , 
for the different instruments. The results obtained for MAGIC (1.3% of Crab flux above 
110 GeV) and CTA (0.30% of Crab flux above 75 GeV) are consistent with those published 
by the respective Collaborations ([15], [24]). On the other hand, our results for VERITAS, 
estimated assuming the A e tf from [27], imply a sensitivity of 0.32% of Crab flux above 165 
GeV, more than a factor 2 better than the one reported at [31]. And given how the DM 
constraints reported in [20] are stronger than those computed in this work, we can expect 
this discrepancy to be even larger for the A e ff actually used in the analysis by the VERITAS 
Collaboration. 

4.3 Annihilation into 77 

Although theoretically highly suppressed, direct annihilation of the DM particles into gamma 
rays would result in a presence of a sharp line in the energy spectrum, whose detection would 
represent an unambiguous confirmation of the DM existence. 

We have made estimates, using the full likelihood approach, of the sensitivity MAGIC 
and CTA observatories have to spectral lines, for the DM particle mass in the energy range 
between 100 GeV and 5 TeV. We require a 5a detection in 50 hours of observations of Segue 1, 
with Einasto[32] DM profile (J = 1.7 x 10 19 GeV 2 cm -5 [33]). For the size of the background 
region we take r = 12. Results are shown on figure 9. MAGIC and CTA exhibit greatest 
sensitivity for m x around 200 and 500 GeV, respectively, with CTA being a factor between 
~5 and ~10 better than MAGIC. Furthermore, for the greater part of the considered m x 
space, the CTA even slightly probes the {av) region below the weak-scale cross-section value, 
~3x 10~ 26 cm 3 s -1 . 

Recent work from Weniger (2012) [21] claims a hint of a gamma-ray line at energy of 
~129 GeV, from 43 months of Galactic Halo public Fermi data and search restricted to the 
20 - 300 GeV range. The inferred value of {av) is 1.27xl0~ 27 cm 3 s -1 (assuming Einasto 
profile). As seen from figure 9, neither MAGIC nor CTA are close to reaching that sensitivity 
using Segue 1 observations. 

On the other hand, observations of a source of much higher J should result in better 
constraints. For example, 50 hours of CTA observations of the Galactic Halo (assuming 
NFW density profile, J = 3.3 x 10 21 GeV 2 cm -5 and r = 2) would yield order of 30 times 
better sensitivity than in the case of Segue 1 (figure 9). Furthermore, 50 hours of data 
would be sufficient for the CTA to test Weniger's claim. It must be noted, however, that this 
computation does not take into account the systematic uncertainties, which might be relevant 
for this search (given how the gamma-ray rate would be ~2% that of the background). 

5 dN/dE = 5.8 x i - 1 3( J B/300GeV)- 2 - 32 - - 131 ° g io (E/300GcV) GeV" 1 cm- 2 s" 1 [15] 
6 Defined as the peak of gamma-ray rate true energy distribution 
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Figure 9: MAGIC (blue) and CTA (red) sensitivities on the (erv) values as a function of m x , for 50 hours of 
Segue 1 observations, assuming a 100% branching ratio into 77. Also plotted are the sensitivity predictions 
for the CTA observations of the Galactic Halo (orange), as well as the (av) value corresponding to DM signal 
hint, at m x — 129 GeV, claimed in [21]. 



5 Discussion 

We have presented an analysis approach for IACTs that uses the full likelihood method, 
constructed to take the maximal advantage of the unique spectral features of DM origin. 
Almost solely through the inclusion of the a priori knowledge on the expected gamma-ray 
spectrum in the likelihood, this method accedes better sensitivity of the analysis, with Im- 
provement Factors reaching values up to order of 10 (depending on the signal model) with 
respect to the recent IACT results. In addition, as shown in section 3.4, these improvements 
are rather insensitive to other analysis characteristics, like the background estimation or 
signal-to-background discrimination. As a result, the full likelihood can be combined to any 
other analysis developments aimed at further sensitivity enhancements. 

In this work, we have focused on the indirect searches for DM annihilation signals with 
IACTs. This is reflected in the specific form of the likelihood function (eq.(2.2)), determined 
by the fact that IACT observations are pointed, cover a relatively narrow field of view, and 
are dominated by background events. Although, to our knowledge, never used for IACTs, 
this concept is a well known analysis method, successfully applied in other fields, including 
DM searches with different techniques and instruments. For instance, a similar approach is 
employed in the direct detection experiments, like XENON100 [34], and even more exten- 
sively, in the indirect searches for DM signals in gamma rays by the Fermi-LAT 7 (see, e.g., 
[35-37]). 

The proposed method is sufficiently general to be used in studies of other physics cases 
studied by the IACTs, the only condition being that a prediction about the expected spectral 

7 http://fermi. gsfc.nasa.gov 
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distribution can be made. A trivial example is the search for DM decay signals, for which we 
only need to substitute the p 2 term by p in eq.(4.3) and (crv) /2m x by l/r x m x in eq.(4.2) (with 
t x referring to the DM particle decay lifetime). Another example where the full likelihood 
can be successfully applied is in the search of the AGN spectra for signatures induced by 
the oscillations of gammas into axion-like particles in the presence of intergalactic magnetic 
fields [38]. This case, however, would require the a priori assumptions on the AGN emission 
and effects of gamma rays interacting with the extragalactic background light. 

Very important characteristic of the full likelihood method (and any likelihood function- 
based analysis) is that it allows a rather straightforward combination of the results obtained 
by different instruments and from different targets. For a given DM model M(6), and A^; nst 
different instruments (or measurements), a global likelihood function can be simply written 
as: 

C T (M{6)) = H d(M(0)). (5.1) 

i=l 

This approach eliminates the complexity required for a common treatment of data and re- 
sponse functions of different telescopes or analyses, required by, e.g. the data stacking method 
(see, e.g. [39]). Within the likelihood scheme, the details of each experiment do not need 
to be combined or averaged. The only necessary information is the value of the likelihood, 
expressed as a function of the free parameter (e.g. {crv)) of a given model for different in- 
struments. Since DM signals are universal and do not depend on the observed target, the 
results from different sources can also be combined through the overall likelihood function 
(as done by Fermi [40]), providing therefore a more sensitive DM search. For example, com- 
bined results (of similar sensitivities to (crv)) from three different observatories (e.g. MAGIC, 
VERITAS and HESS 8 ) would benefit from an extra improvement in the sensitivity by a fac- 
tor of ~ 70%. This approach would offer the best chances of discovering DM in indirect 
VHE gamma-ray searches or of setting the most stringent limits attainable by this kind of 
observations, placing therefore a new landmark in the field. 
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